Este estudio presenta un análisis básico de RNA-Seq con énfasis en las posibilidades de visualización de datos y resultados^.
En su versión actual los gráficos son correctos pero poco elaborados.
Para mejorarlos se utilizaran gráficos con ggplot2
utilizando como modelo algunas de las propuestas del documento: https://github.com/UofABioinformaticsHub/DataVisualisaton_BIS2016.git
Este estudio se llevará a cabo usando R/Bioconductor por lo que es preciso tener instaladas un conjunto de librerías. Esto puede hacerse siguiendo el procedimiento descrito a continuación.
El código que se presenta debería ejecutarse tan sólo una vez!
if (!require(BiocManager)) {
install.packages("BiocManager", dep = TRUE)
}
installifnot <- function(pckgName, BioC = TRUE) {
if (BioC) {
if (!require(pckgName, character.only = TRUE)) {
BiocManager::install(pckgName)
}
} else {
if (!require(pckgName, character.only = TRUE)) {
install.packages(pckgName, dep = TRUE)
}
}
}
installifnot("limma")
installifnot("edgeR")
installifnot("org.Hs.eg.db")
installifnot("clusterProfiler")
installifnot("dplyr", BioC = FALSE)
installifnot("gplots", BioC = FALSE)
installifnot("ggvenn", BioC = FALSE)
installifnot("pheatmap")
installifnot("prettydoc", BioC = FALSE)
installifnot("ggnewscale", BioC = FALSE)
# Pels grafics que empren ggplot2
installifnot("locfit", BioC = FALSE)
installifnot("magrittr", BioC = FALSE)
installifnot("statmod", BioC = FALSE)
# Pel HeatMap Interactiu
installifnot("heatmaply")
# Pel VolcanoPlot interactiu
installifnot("here")
installifnot("janitor") # Cleaning column names
installifnot("scales") # Transform axis scales
installifnot("ggrepel")
Puede resultar útil, aunque no es para nada imprescindible, disponer
de, al menos dos, directorios específicos: - datos (o un
nombre similar) donde guardar y de donde cargar los datos. -
resultso (o un nombre similar) donde escribir los
resultados.
if (!dir.exists("datos")) dir.create("datos")
if (!dir.exists("results")) dir.create("results")
if (!dir.exists("figures")) dir.create("figures")
Este análisis se basa en un estudio publicado recientemente (Arunachalam 2020) que investigaba la respuesta inmune a la infección con SARS-COV-2 desde un a perspectiva de biología de sistemas utilizando tecnología de secuenciación, RNA-Seq.
Los datos se han depositado en el repositorio “Gene Expression Omnibus (GEO)” con el identificador GSE152418.
Los datos de contajes se han descargado desde l repositorio “GEO”
al archivo Rawcounts.csv que se utilizará como punto de
partida para el análisis.
counts <- read.csv("datos/RawCounts.csv", row.names = 1)
selectedCounts <- as.matrix(counts)
En la misma web se dispone de información los grupos y otras covariables o variables auxiliares. Con éstos se crea un objeto (habitualmente denominado “targets”) que debe estar sincronizado con el anterior, es decir, cuyas filas deben corresponderse con las columnas de la matriz de datos.
sampleNames <- colnames(selectedCounts)
grupos <- c(rep("COVID", 17), rep("SANO", 17))
colores = c(rep("red", 17), rep("blue", 17))
selectedTargets <- data.frame(samples = sampleNames, group = grupos, cols = colores)
rownames(selectedTargets) <- selectedTargets[, 1]
Los datos de secuenciación pueden estar “desbalanceados” en el sentido que distintas muestras pueden contener un número distinto de secuencias, lo que puede inducir erróneamente a pensar que un gen se expresa más en una muestra que en otra, cuando esto se deba a esta diferencia global.
Esto puede evitarse expresando los contajes como “CPMs” es decir “counts per million”, lo que no modificará los resultados de comparaciones posteriores, pero hará que las muestras sean comparables en número , lo que es útil y necesario para los análisis posteriores.
library(edgeR)
selectedCounts[1:5, 1:6]
COV145 COV147 COV149 COV150 COV151 COV152
ENSG00000223972 0 0 0 0 0 0
ENSG00000227232 1 0 3 16 1 2
ENSG00000278267 3 1 2 0 3 0
ENSG00000243485 2 0 0 1 0 0
ENSG00000284332 0 0 0 0 0 0
apply(selectedCounts, 2, sum)
COV145 COV147 COV149 COV150 COV151 COV152 COV153 COV154
15048983 13946337 13520020 18863644 13659353 13279390 12967503 14300187
COV155 COV156 COV157 COV175 COV176 COV177 COV178 COV179
17192277 15626012 17041245 13415720 13679816 16092084 17379155 13672179
COV180 HEA061 HEA062 HEA063 HEA064 HEA065 HEA066 HEA067
13722892 17003281 17359694 15570080 16172199 19420885 18264376 16398893
HEA068 HEA069 HEA070 HEA071 HEA181 HEA182 HEA183 HEA184
14895860 14939005 20458479 17820834 14973734 14280768 18297901 17177433
HEA185 HEA186
16467873 16860317
counts.CPM <- cpm(selectedCounts)
counts.CPM[1:5, 1:6]
COV145 COV147 COV149 COV150 COV151 COV152
ENSG00000223972 0.00000000 0.00000000 0.0000000 0.00000000 0.00000000 0.0000000
ENSG00000227232 0.06644967 0.00000000 0.2218932 0.84819243 0.07320991 0.1506093
ENSG00000278267 0.19934902 0.07170342 0.1479288 0.00000000 0.21962973 0.0000000
ENSG00000243485 0.13289935 0.00000000 0.0000000 0.05301203 0.00000000 0.0000000
ENSG00000284332 0.00000000 0.00000000 0.0000000 0.00000000 0.00000000 0.0000000
apply(counts.CPM, 2, sum)
COV145 COV147 COV149 COV150 COV151 COV152 COV153 COV154 COV155 COV156 COV157
1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06
COV175 COV176 COV177 COV178 COV179 COV180 HEA061 HEA062 HEA063 HEA064 HEA065
1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06
HEA066 HEA067 HEA068 HEA069 HEA070 HEA071 HEA181 HEA182 HEA183 HEA184 HEA185
1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06
HEA186
1e+06
Una vez los datos estan como CPMs, se procede a filtrarlos,
Los genes con recuentos muy bajos en todas las librerías (es decir en todas las muestras) proporcionan poca evidencia de expresión diferencial, por lo que es habitual eliminar aquellos genes que, o bien varían muy poco entre grupos, o bien presentan poca o nula expresión en la mayoría de las muestras.
En este caso, siguiendo un criterio habitual, se opta por conservar únicamente aquellos genes que presentan algún valor en, al menos, tres muestras de cada grupo.
thresh <- counts.CPM > 0
keep <- (rowSums(thresh[, 1:10]) >= 3) & (rowSums(thresh[, 11:20]) >= 3)
counts.keep <- counts.CPM[keep, ]
dim(counts.CPM)
[1] 60683 34
dim(counts.keep)
[1] 33039 34
Aunque no sea más que un ejemplo basta con ver los dos primeros genes para comprobar como el primero no cumple la condición, en el grupo “SANO”, mientras que los siguientes sí que la cumplen. Por lo tanto, al filtrar desaparece el primer gen de la matriz filtrada, pero no los dos siguientes.
round(head(counts.CPM), 3)
COV145 COV147 COV149 COV150 COV151 COV152 COV153 COV154 COV155
ENSG00000223972 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.00 0.116
ENSG00000227232 0.066 0.000 0.222 0.848 0.073 0.151 1.388 0.42 1.163
COV156 COV157 COV175 COV176 COV177 COV178 COV179 COV180 HEA061
ENSG00000223972 0.000 0.117 0.000 0.000 0.000 0.058 0.000 0.000 0.000
ENSG00000227232 0.000 0.117 0.000 0.950 1.989 0.173 0.146 0.656 0.353
HEA062 HEA063 HEA064 HEA065 HEA066 HEA067 HEA068 HEA069 HEA070
ENSG00000223972 0.000 0.000 0.062 0.000 0.000 0.000 0.000 0.000 0.000
ENSG00000227232 0.000 0.000 0.062 0.103 0.055 0.000 0.067 0.067 0.049
HEA071 HEA181 HEA182 HEA183 HEA184 HEA185 HEA186
ENSG00000223972 0.000 0.0 0.00 0.000 0.000 0.000 0.000
ENSG00000227232 0.112 0.2 0.28 0.055 0.000 0.000 0.415
[ reached getOption("max.print") -- omitted 4 rows ]
round(head(counts.keep), 3)
COV145 COV147 COV149 COV150 COV151 COV152 COV153 COV154 COV155
ENSG00000227232 0.066 0.000 0.222 0.848 0.073 0.151 1.388 0.420 1.163
ENSG00000278267 0.199 0.072 0.148 0.000 0.220 0.000 0.000 0.070 0.174
COV156 COV157 COV175 COV176 COV177 COV178 COV179 COV180 HEA061
ENSG00000227232 0.000 0.117 0.000 0.950 1.989 0.173 0.146 0.656 0.353
ENSG00000278267 0.128 0.117 0.224 0.219 0.435 0.000 0.219 0.000 0.235
HEA062 HEA063 HEA064 HEA065 HEA066 HEA067 HEA068 HEA069 HEA070
ENSG00000227232 0.000 0.000 0.062 0.103 0.055 0.000 0.067 0.067 0.049
ENSG00000278267 0.288 0.193 0.495 0.154 0.493 0.427 0.134 0.469 0.098
HEA071 HEA181 HEA182 HEA183 HEA184 HEA185 HEA186
ENSG00000227232 0.112 0.200 0.28 0.055 0.000 0.000 0.415
ENSG00000278267 0.337 0.200 0.14 0.328 0.699 0.304 0.178
[ reached getOption("max.print") -- omitted 4 rows ]
Cuando se trabaja con distintos objetos referidos a unos mismos datos, como la matriz de contajes y el objeto “targets”, es útil disponer de contenedores que permitan trabajar con todos ellos a la vez, lo que no sólo facilita el trabajo sino que ayuda a evitar “desincronizaciones”.
Éste es el caso de la clase ExpressionSet habitualmente
utilizada con microarrays o de la clase que la generaliza, llamada
SummarizedExperiment.
Para datos de contaje es habitual usar una clase similar a
ExpressionSet llamada DGEList” pensadas para
manejar datos de contajes , definida en el paquete edgeR.
Esta clase, más simple que las anteriores, utiliza listas para almacenar
recuentos de “reads” e información asociada de tecnologías de
secuenciación. Puede encontrarse información al respecto en la ayuda del
paquete edgeR.
dgeObj <- DGEList(counts = counts.keep, lib.size = colSums(counts.keep), norm.factors = rep(1,
ncol(counts.keep)), samples = selectedTargets, group = selectedTargets$group,
genes = rownames(counts.keep), remove.zeros = FALSE)
show(dgeObj)
An object of class "DGEList"
$counts
COV145 COV147 COV149 COV150 COV151
ENSG00000227232 0.06644967 0.00000000 0.2218932 0.84819243 0.07320991
ENSG00000278267 0.19934902 0.07170342 0.1479288 0.00000000 0.21962973
COV152 COV153 COV154 COV155 COV156 COV157
ENSG00000227232 0.15060933 1.3880853 0.41957493 1.163313 0.0000000 0.11736232
ENSG00000278267 0.00000000 0.0000000 0.06992916 0.174497 0.1279917 0.11736232
COV175 COV176 COV177 COV178 COV179 COV180
ENSG00000227232 0.0000000 0.9503052 1.9885554 0.1726206 0.14628246 0.6558384
ENSG00000278267 0.2236183 0.2193012 0.4349965 0.0000000 0.21942369 0.0000000
HEA061 HEA062 HEA063 HEA064 HEA065 HEA066
ENSG00000227232 0.3528731 0.0000000 0.0000000 0.06183451 0.10298192 0.05475139
ENSG00000278267 0.2352487 0.2880235 0.1926772 0.49467608 0.15447288 0.49276252
HEA067 HEA068 HEA069 HEA070 HEA071 HEA181
ENSG00000227232 0.0000000 0.06713275 0.06693886 0.04887949 0.1122282 0.2003508
ENSG00000278267 0.4268581 0.13426549 0.46857204 0.09775898 0.3366846 0.2003508
HEA182 HEA183 HEA184 HEA185 HEA186
ENSG00000227232 0.28009698 0.05465108 0.0000000 0.0000000 0.41517606
ENSG00000278267 0.14004849 0.32790646 0.6985910 0.3036215 0.17793260
[ reached getOption("max.print") -- omitted 3 rows ]
33034 more rows ...
$samples
group lib.size norm.factors samples group.1 cols
COV145 COVID 999825.8 1 COV145 COVID red
COV147 COVID 999530.1 1 COV147 COVID red
COV149 COVID 999667.8 1 COV149 COVID red
COV150 COVID 999872.5 1 COV150 COVID red
COV151 COVID 999904.2 1 COV151 COVID red
29 more rows ...
$genes
genes
ENSG00000227232 ENSG00000227232
ENSG00000278267 ENSG00000278267
ENSG00000233750 ENSG00000233750
ENSG00000268903 ENSG00000268903
ENSG00000269981 ENSG00000269981
33034 more rows ...
Uno de los aspectos interesantes de estas clases es la posibilidad de extraer partes de todos los objetos a la vez con el operador de “subsetting”.
dim(dgeObj)
[1] 33039 34
dgeObjShort <- dgeObj[, c(1:5, 11:15)]
# Library size information is stored in the samples slot
dgeObjShort$samples
group lib.size norm.factors samples group.1 cols
COV145 COVID 999825.8 1 COV145 COVID red
COV147 COVID 999530.1 1 COV147 COVID red
COV149 COVID 999667.8 1 COV149 COVID red
COV150 COVID 999872.5 1 COV150 COVID red
COV151 COVID 999904.2 1 COV151 COVID red
COV157 COVID 999582.4 1 COV157 COVID red
COV175 COVID 998800.7 1 COV175 COVID red
COV176 COVID 999667.6 1 COV176 COVID red
COV177 COVID 999807.4 1 COV177 COVID red
COV178 COVID 999721.9 1 COV178 COVID red
colnames(dgeObjShort$counts)
[1] "COV145" "COV147" "COV149" "COV150" "COV151" "COV157" "COV175" "COV176"
[9] "COV177" "COV178"
# https://hbctraining.github.io/DGE_workshop_salmon/lessons/AnnotationDbi_lesson.html
annotations_orgDb <- AnnotationDbi::select(org.Hs.eg.db, # database
keys = dgeObj$genes[,1], # data to use for retrieval
columns = c("SYMBOL", "ENTREZID","GENENAME"), # information to retreive for given data
keytype = "ENSEMBL")
length(which(is.na(annotations_orgDb$SYMBOL) == FALSE))
[1] 23304
non_duplicates_idx <- which(duplicated(annotations_orgDb$SYMBOL == FALSE))
annotations_orgDb<- annotations_orgDb[non_duplicates_idx, ]
# https://hbctraining.github.io/Intro-to-R-flipped/lessons/09_reordering-to-match-datasets.html
cbind(dgeObj$genes[, 1], annotations_orgDb$ENSEMBL)[1:100, ]
[,1] [,2]
[1,] "ENSG00000227232" "ENSG00000278267"
[2,] "ENSG00000278267" "ENSG00000233750"
[3,] "ENSG00000233750" "ENSG00000269981"
[4,] "ENSG00000268903" "ENSG00000241860"
[5,] "ENSG00000269981" "ENSG00000279928"
[6,] "ENSG00000241860" "ENSG00000279457"
[7,] "ENSG00000279928" "ENSG00000228463"
[8,] "ENSG00000279457" "ENSG00000237094"
[9,] "ENSG00000228463" "ENSG00000230021"
[10,] "ENSG00000237094" "ENSG00000225972"
[11,] "ENSG00000230021" "ENSG00000225630"
[12,] "ENSG00000225972" "ENSG00000237973"
[13,] "ENSG00000225630" "ENSG00000229344"
[14,] "ENSG00000237973" "ENSG00000248527"
[15,] "ENSG00000229344" "ENSG00000198744"
[16,] "ENSG00000248527" "ENSG00000228327"
[17,] "ENSG00000198744" "ENSG00000237491"
[18,] "ENSG00000228327" "ENSG00000230092"
[19,] "ENSG00000237491" "ENSG00000177757"
[20,] "ENSG00000230092" "ENSG00000228794"
[21,] "ENSG00000177757" "ENSG00000225880"
[22,] "ENSG00000228794" "ENSG00000288531"
[23,] "ENSG00000225880" "ENSG00000234711"
[24,] "ENSG00000288531" "ENSG00000272438"
[25,] "ENSG00000234711" "ENSG00000230699"
[26,] "ENSG00000272438" "ENSG00000223764"
[27,] "ENSG00000230699" "ENSG00000187634"
[28,] "ENSG00000223764" "ENSG00000188976"
[29,] "ENSG00000187634" "ENSG00000187961"
[30,] "ENSG00000188976" "ENSG00000187583"
[31,] "ENSG00000187961" "ENSG00000187642"
[32,] "ENSG00000187583" "ENSG00000272512"
[33,] "ENSG00000187642" "ENSG00000188290"
[34,] "ENSG00000272512" "ENSG00000187608"
[35,] "ENSG00000188290" "ENSG00000224969"
[36,] "ENSG00000187608" "ENSG00000188157"
[37,] "ENSG00000224969" "ENSG00000217801"
[ reached getOption("max.print") -- omitted 63 rows ]
reorder_idx <- match(dgeObj$genes[, 1], annotations_orgDb$ENSEMBL)
annotations_orgDb_reordered <- annotations_orgDb[reorder_idx, ]
cbind(dgeObj$genes[, 1], annotations_orgDb_reordered$ENSEMBL)[1:100, ]
[,1] [,2]
[1,] "ENSG00000227232" NA
[2,] "ENSG00000278267" "ENSG00000278267"
[3,] "ENSG00000233750" "ENSG00000233750"
[4,] "ENSG00000268903" NA
[5,] "ENSG00000269981" "ENSG00000269981"
[6,] "ENSG00000241860" "ENSG00000241860"
[7,] "ENSG00000279928" "ENSG00000279928"
[8,] "ENSG00000279457" "ENSG00000279457"
[9,] "ENSG00000228463" "ENSG00000228463"
[10,] "ENSG00000237094" "ENSG00000237094"
[11,] "ENSG00000230021" "ENSG00000230021"
[12,] "ENSG00000225972" "ENSG00000225972"
[13,] "ENSG00000225630" "ENSG00000225630"
[14,] "ENSG00000237973" "ENSG00000237973"
[15,] "ENSG00000229344" "ENSG00000229344"
[16,] "ENSG00000248527" "ENSG00000248527"
[17,] "ENSG00000198744" "ENSG00000198744"
[18,] "ENSG00000228327" "ENSG00000228327"
[19,] "ENSG00000237491" "ENSG00000237491"
[20,] "ENSG00000230092" "ENSG00000230092"
[21,] "ENSG00000177757" "ENSG00000177757"
[22,] "ENSG00000228794" "ENSG00000228794"
[23,] "ENSG00000225880" "ENSG00000225880"
[24,] "ENSG00000288531" "ENSG00000288531"
[25,] "ENSG00000234711" "ENSG00000234711"
[26,] "ENSG00000272438" "ENSG00000272438"
[27,] "ENSG00000230699" "ENSG00000230699"
[28,] "ENSG00000223764" "ENSG00000223764"
[29,] "ENSG00000187634" "ENSG00000187634"
[30,] "ENSG00000188976" "ENSG00000188976"
[31,] "ENSG00000187961" "ENSG00000187961"
[32,] "ENSG00000187583" "ENSG00000187583"
[33,] "ENSG00000187642" "ENSG00000187642"
[34,] "ENSG00000272512" "ENSG00000272512"
[35,] "ENSG00000188290" "ENSG00000188290"
[36,] "ENSG00000187608" "ENSG00000187608"
[37,] "ENSG00000224969" "ENSG00000224969"
[ reached getOption("max.print") -- omitted 63 rows ]
dgeObj$genes <- annotations_orgDb_reordered
dgeObj$genes[1:100, ]
ENSEMBL SYMBOL ENTREZID
NA <NA> <NA> <NA>
2 ENSG00000278267 MIR6859-1 102466751
3 ENSG00000233750 CICP27 100420257
NA.1 <NA> <NA> <NA>
5 ENSG00000269981 <NA> <NA>
6 ENSG00000241860 <NA> <NA>
7 ENSG00000279928 <NA> <NA>
8 ENSG00000279457 WASH9P 102723897
9 ENSG00000228463 RPL23AP21 728481
10 ENSG00000237094 <NA> <NA>
11 ENSG00000230021 LOC101928626 101928626
12 ENSG00000225972 MTND1P23 100887749
13 ENSG00000225630 MTND2P28 100652939
14 ENSG00000237973 MTCO1P12 107075141
15 ENSG00000229344 MTCO2P12 107075310
16 ENSG00000248527 MTATP6P1 106480796
17 ENSG00000198744 MTCO3P12 107075270
18 ENSG00000228327 <NA> <NA>
GENENAME
NA <NA>
2 microRNA 6859-1
3 capicua transcriptional repressor pseudogene 27
NA.1 <NA>
5 <NA>
6 <NA>
7 <NA>
8 WAS protein family homolog 9, pseudogene
9 ribosomal protein L23a pseudogene 21
10 <NA>
11 uncharacterized LOC101928626
12 MT-ND1 pseudogene 23
13 MT-ND2 pseudogene 28
14 MT-CO1 pseudogene 12
15 MT-CO2 pseudogene 12
16 MT-ATP6 pseudogene 1
17 MT-CO3 pseudogene 12
18 <NA>
[ reached 'max' / getOption("max.print") -- omitted 82 rows ]
Aunque podríamos haber creado el objeto a partir de todas las muestras, y haber realizado la extracción de genes y muestras posteriormente, hemos optado por no hacerlo para facilitar el seguimiento del proceso.
Además de estandarizar los contajes, es importante eliminar otros sesgos de composición entre librerías. Esto puede hacerse aplicando la normalización por el método TMM que genera un conjunto de factores de normalización, tal que producto de estos factores y los tamaños de librería (el número de secuencias de cada muestra) definen el tamaño efectivo de dichas muestras, es decir el peso real que se les asignará en las comparaciones posteriores.
Aunque esto puede parecer artificial, no lo es porque la normalización tiene en cuenta otros factores, como el sesgo de composición entre librerías, que podrían hacer que los mismos valores en distintas muestras no reflejaran su importancia relativa.
La función calcNormFactors, de la librería
edgeR, calcula los factores de normalización
mencionados.
library(edgeR)
dgeObj_norm <- calcNormFactors(dgeObj)
head(dgeObj_norm$samples, 10)
group lib.size norm.factors samples group.1 cols
COV145 COVID 999825.8 1.0018278 COV145 COVID red
COV147 COVID 999530.1 1.1010464 COV147 COVID red
COV149 COVID 999667.8 0.9553476 COV149 COVID red
COV150 COVID 999872.5 0.9479175 COV150 COVID red
COV151 COVID 999904.2 0.8996072 COV151 COVID red
COV152 COVID 999763.7 1.0033710 COV152 COVID red
COV153 COVID 999874.5 0.9870834 COV153 COVID red
COV154 COVID 999736.9 0.8538699 COV154 COVID red
COV155 COVID 999399.0 0.8230210 COV155 COVID red
COV156 COVID 999868.7 0.8878678 COV156 COVID red
Esto no modificará la matriz de contajes, pero actualizará los factores de normalización en el objeto DGEList (sus valores predeterminados son 1).
head(dgeObj$samples)
group lib.size norm.factors samples group.1 cols
COV145 COVID 999825.8 1 COV145 COVID red
COV147 COVID 999530.1 1 COV147 COVID red
COV149 COVID 999667.8 1 COV149 COVID red
COV150 COVID 999872.5 1 COV150 COVID red
COV151 COVID 999904.2 1 COV151 COVID red
COV152 COVID 999763.7 1 COV152 COVID red
head(dgeObj_norm$samples)
group lib.size norm.factors samples group.1 cols
COV145 COVID 999825.8 1.0018278 COV145 COVID red
COV147 COVID 999530.1 1.1010464 COV147 COVID red
COV149 COVID 999667.8 0.9553476 COV149 COVID red
COV150 COVID 999872.5 0.9479175 COV150 COVID red
COV151 COVID 999904.2 0.8996072 COV151 COVID red
COV152 COVID 999763.7 1.0033710 COV152 COVID red
Es decir, aunque no se observen cambios en la matriz de contajes, cuando se utilizan estos factores de normalización en algún cálculo la importancia de las distintas columnas se tendrá en cuenta.
Las transformaciones anteriores buscan compensar el tamaño distinto de las librerías o la distinta composición de éstas, pero las distribuciones de los contajes en cada muestra son asimétricas.
boxplot(dgeObj_norm$counts, col = dgeObj_norm$samples$cols, las = 2, cex.axis = 0.7,
main = "Contajes normalizados", ylim = c(0, 10000))
Para finalizar el preprocesado se toman logaritmo base dos de los contajes.
log2count_norm <- cpm(dgeObj_norm, log = TRUE)
source("Rcode/niceBoxPlot.R")
gg1 <- niceBoxPlot(log2count_norm)
ggsave(filename = "figures/boxplot.png", gg1, bg = "transparent") # ,dpi = 300
gg1
También se puede hacer usando plotly.
niceBoxPlotly(log2count_norm)